Mini pruebas:

1 ciclo de asimilación dentro del experimento 20181122, para el 20181123000000 donde hay disponible observaciones del AMSU-A en NOAA15 y NOAA 18.

files <- list.files("analisis/radianzas", pattern = "amsua_n15", full.names = TRUE) %>%
  .[!str_detect(., "conv")]
diag <- map(files, function(f){
  meta <- unglue::unglue(f, "analisis/radianzas/asim_{sensor}_{plat}_{date}.ensmean_{exp}")
  # print(f)
  out <- fread(f)
  # .[V10 == 1] %>% 
  
  if (file.size(f) != 0) {
    out[, date := ymd_hms(meta[[1]][["date"]])] %>% 
      .[, exp := meta[[1]][["exp"]]]
  }
  out
}) %>%
  rbindlist()

colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
                    "errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld", 
                    "rcldp", paste0("pred", seq(8)), "date", "exp")
diag[qc == 0, .N, by = .(channel, exp)] %>% 
  ggplot(aes(factor(channel), exp)) +
  geom_raster(aes(fill = N)) +
  scale_fill_viridis_c() +
  labs(x = "channel", y = "", subtitle = "QC = 0") +
  theme_minimal()

diag[qc == 0 & errinv !=0, .N, by = .(channel, exp)] %>% 
  dcast(exp ~ channel) %>% 
  knitr::kable(caption = "Cantidad de observaciones potencialmente asimilables (QC == = y errinv != 0) para cada prueba y cada canal del AMSU-A")
## Using 'N' as value column. Use 'value.var' to override
Cantidad de observaciones potencialmente asimilables (QC == = y errinv != 0) para cada prueba y cada canal del AMSU-A
exp 1 2 3 4 5 6 7 8 15
QC4new 127 129 129 129 129 404 759 420 126
QC4old 105 107 107 107 107 404 759 421 104
QC4ori 103 105 105 105 105 404 759 420 102
QC6ecmwf 99 101 101 101 101 403 709 410 98
QC6ecmwfall 184 183 176 186 192 402 716 413 181
allsky 147 146 139 149 155 328 689 396 144
clearsky 65 67 67 67 67 328 679 389 64
noQC6 99 101 101 101 101 496 709 410 98 
diag[channel %in% c(1:8)] %>% 
  dcast(lon + lat + channel ~ exp, value.var = "qc") %>% 
  .[, diff := allsky - clearsky] %>% 
  .[diff != 0] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(color = "darkorange", size = 0.7) +
  # geom_point(aes(color = errinv)) +
  scale_color_viridis_c() +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_grid(.~ channel) +
  labs(x = "lon",
       subtitle = "Ubicación de observaciones extras en allsky respecto de clearsky") +
  theme_minimal()

diag[channel %in% c(4)] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = factor(abs(qc))), size = 0.7) +
  scale_color_viridis_d() +
   geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  # facet_grid(exp~channel) +
  facet_wrap(~exp, ncol = 4) +
  labs(subtitle = "Control de calidad para el canal 4",
       x = "", y = "",
       color = "QC") +
  theme_minimal()

diag[channel %in% c(4:8) & qc == 0 & errinv != 0] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = tbc), size = 0.7) +
  scale_color_viridis_c(option = "A") +
  # scale_color_paletteer_c("ggthemes::Red-Blue Diverging") +
   geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_grid(channel~exp) +
  labs(x = "", y = "",
       color = "O-B") +
  theme_minimal()

# diag[channel %in% c(1:8)] %>% 
#   ggplot(aes(ConvertLongitude(lon), lat)) +
#   geom_point(aes(color = factor(qc))) +
#   scale_color_viridis_d() +
#    geom_mapa() +
#   coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
#   facet_grid(exp~channel) +
#   theme_minimal()


# diag[channel %in% c(1:15) & errinv != 0, .N, by = .(channel, exp)] %>% 
#   dcast(exp~channel)

diag[, qc := abs(qc)] %>% 
  .[channel %in% (5), .N, by = .(exp, qc)] %>% 
  dcast(exp~qc) %>% 
  knitr::kable(caption = "Categorías de control de calidad por experimento para el canal 5.")
## Using 'N' as value column. Use 'value.var' to override
Categorías de control de calidad por experimento para el canal 5.
exp 0 3 4 8 50 51
QC4new 129 NA 25 228 308 92
QC4old 107 NA 23 138 307 207
QC4ori 105 NA 21 139 308 209
QC6ecmwf 101 NA 18 142 310 211
QC6ecmwfall 192 3 21 87 311 168
allsky 155 3 4 38 450 132
clearsky 67 NA 3 91 450 171
noQC6 101 NA 34 190 NA 457 
diag[channel %in% c(1:15) & qc == 0 & errinv > 0.000031623,
     .N, by = .(exp)] %>% 
  dcast(exp~.) %>% 
  .[order(-.)] %>% 
  knitr::kable(caption = "Cantidad total de observaciones potencialmente asimilables")
## Using 'N' as value column. Use 'value.var' to override
Cantidad total de observaciones potencialmente asimilables
exp .
QC6ecmwfall 2633
QC4new 2352
allsky 2293
QC4old 2221
noQC6 2216
QC4ori 2208
QC6ecmwf 2123
clearsky 1793 
diag[channel %in% c(4:8)] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = errinv)) +
  scale_color_viridis_c() +
   geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_grid(channel~exp) +
  labs(x = "", y = "") +
  theme_minimal()

files <- list.files("analisis/radianzas", pattern = "analysis", full.names = TRUE) %>% 
    .[!str_detect(., "asim")]

update <- map(files, function(f) {

  descriptores <- unglue::unglue(f, c("analisis/radianzas/analysis.ensmean_{exp}"))

  ana <- ReadNetCDF(f, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                 lon = "XLONG", lat = "XLAT")) %>%
    .[, p := p + PB] %>%
    .[, t := tk(t, p)] %>%
    .[, rh := rh(qv, p, t)] %>% 
    .[, td := td(qv, p) + 273.15] %>% 
    .[, ":="(Time = NULL,
             # west_east = NULL,
             # south_north = NULL,
             qv = NULL,
             PB = NULL)] %>% 
    .[, c("u", "v") := uvmet(f)] %>% 
    .[, ":="(exp = descriptores[[1]][["exp"]])] %>% 
    .[]

}) %>% 
  rbindlist()

  guess <- ReadNetCDF("analisis/radianzas/wrfarw.ensmean", 
                      vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                 lon = "XLONG", lat = "XLAT")) %>%
    .[, p := p + PB] %>%
    .[, t := tk(t, p)] %>%
    .[, rh := rh(qv, p, t)] %>% 
    .[, td := td(qv, p) + 273.15] %>% 
    .[, ":="(Time = NULL,
             # west_east = NULL,
             # south_north = NULL,
             qv = NULL,
             PB = NULL)] %>% 
    .[, c("u", "v") := uvmet("analisis/radianzas/wrfarw.ensmean")] %>% 
    # .[, ":="(exp = "guess")] %>% 
    .[]

  
update <- guess[update, on = c("lon", "lat", "bottom_top", "south_north", "west_east")]

Comparación de update para QC4

update[, .(mean_ana = mean(i.t, na.rm = TRUE),
              mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, west_east, exp)] %>% 
  .[exp %in% c("QC4new", "QC4old")] %>% 
  ggplot(aes(west_east, bottom_top)) +
  geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
  scale_fill_divergent() +
  labs(fill = "update T") +
  facet_wrap(~ exp) +
  theme_minimal()

update[, .(mean_ana = mean(i.t, na.rm = TRUE),
              mean_guess = mean(t, na.rm = TRUE)), by = .(bottom_top, south_north, exp)] %>% 
    .[exp %in% c("QC4new", "QC4old")] %>% 
  ggplot(aes(south_north, bottom_top)) +
  geom_contour_fill(aes(z = (mean_ana - mean_guess))) +
  scale_fill_divergent() +
  labs(fill = "update T") +
  facet_wrap(~ exp) +
  theme_minimal()

La diferencia entre ambas pruebas para la temperatura es en promedio del ordende 10-3 o 10-4 y un poco menor en los niveles mas altos. El rango de la diferencia está aproximadamente entre -0.4 y 0.2 K.

update %>% 
  dcast(bottom_top + south_north + west_east ~ exp, value.var = "i.t") %>% 
  .[, .(diff = mean(abs(QC4new - QC4old)),
        max = max(QC4new - QC4old),
        min = min(QC4new - QC4old)), by = .(bottom_top)] %>% 
  melt(measure.vars = c("diff", "max", "min")) %>% 
  ggplot(aes(value, bottom_top)) +
  geom_vline(xintercept = 0, color = "darkgray") +
    geom_path(aes(color = variable)) +
    scale_color_viridis_d() +
  labs(y = "sigma level", color = "") +
  theme_minimal()

Analisis de oservaciones potencialmente asimilables en todo el experimento

files <- list.files("analisis/diagfiles/E7", pattern = "asim", full.names = TRUE) %>%
  .[!str_detect(., "conv")]
diag <- map(files, function(f){
  meta <- unglue::unglue(f, "analisis/diagfiles/{exp}/asim_{sensor}_{plat}_{date}.ensmean")
  # print(f)
  out <- fread(f)
  # .[V10 == 1] %>% 
  
  if (file.size(f) != 0) {
    out[, date := ymd_hms(meta[[1]][["date"]])]
  }
  out
}) %>%
  rbindlist()

colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
                    "errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld", 
                    "rcldp", paste0("pred", seq(8)), "date")

satinfo <- fread("global_satinfo.txt") %>% 
  setnames(old = c("!sensor/instr/sat", "chan"), new = c("sensor", "channel"))
diag[, .N, by = .(sensor, channel, qc)] %>% 
  satinfo[., on = c("sensor", "channel")] %>% 
  ggplot(aes(factor(channel), factor(abs(qc)))) +
  geom_tile(aes(fill = N)) +
  scale_fill_viridis_c() +
  geom_point(data = ~.x[iuse != 1], shape = 4, color = "grey10") +
  # scale_shape_manual(values = c(0, 4)) +
  facet_wrap(~sensor, scales = "free_x") +
  labs(x = "channel", y = "QC",
       caption = "Con 'x' se marcan los canales rechazados por satinfo") +
  theme_minimal()